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^^ I Abstract: We investigate the effects of the saturation boundary on smaU-x evolution at 

the next-to-leading order accuracy and beyond. We demonstrate that the instabilities of 

f~^ ' the next-to-leading order BFKL evolution are not cured by the presence of the nonlinear 

saturation effects, and a resummation of the higher order corrections is therefore needed 
for the nonlinear evolution. The renormalization group improved resummed equation in 
the presence of the saturation boundary is investigated, and the corresponding saturation 
scale is extracted. A significant reduction of the saturation scale is found, and we observe 
C^ \ that the onset of the saturation corrections is delayed to higher rapidities. This seems to be 

related to the characteristic feature of the resummed splitting function which at moderately 
small values of x possesses a minimum. 
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1. Introduction 

It is well known that the BFKL evolution [1-3] suffers, beyond the leading logarithmic 
order (LL), from large corrections which are related to the running of the coupling and 
to kinematical effects, such as energy-momentum conservation. Such corrections can be 
easily added to the leading order formalism in phenomenological applications but it is 
also desirable to have a good control of the next-to-leading logarithmic (NLL) corrections 
in the precise theoretical formulation. BFKL describes an evolution in rapidity where 
the asymptotic limit s/t — )• oo (with s the cms energy and t the momentum transfer) 
described by the LL approximation does not correspond to a vanishing coupling strength a^, 
unlike the QCD renormalization group (RG) evolution dictated by the DGLAP equations 
where the asymptotic limit Q^ — ;■ oo indeed implies that as{Q) —5- 0. Consequently, the 
higher order corrections in the latter are well controlled in the region of the applicability 
of the formalism. In the BFKL framework on the other hand, as s becomes very large. 



the evolution receives contributions from an increasing phase space in momenta where 
as can be very large, and as a consequence there is no reason to expect that the higher 
order corrections can be neglected. Indeed, the calculated NLL corrections to the BFKL 
formalism turn out to be very large [4-9], and it is therefore important to try include the 
effects of these large corrections in any application. 

The NLL BFKL evolution is affected by certain problems which lead to unstable results 
such as negative and oscillating cross sections, see for example [10]. These instabilities of 
the NLL evolution originate from the existence of negative double and triple poles in the 
eigenvalue of the evolution kernel. Certain strategies have therefore been proposed to 
deal with the instabilities of the formalism. These involve the all-order resummation of 
the dominant parts of the higher order corrections, and there exist several prescriptions 
[11-20], all of which are consistent with each other. These procedures lead to the so-called 
"renormalization group improved" small-x evolution. For a nice and comprehensive review 
of the NLL formalism, its problems, and the resummation strategies, see [10]. 

When the QCD dynamics is probed at very small-x, however, we also expect corrections 
from so-called saturation effects which are related to the formation of strong classical 
color fields. The effective theory which takes into account these effects is the Color Glass 
Condensate (CGC) (for a review see [21]). These corrections lead to a generalization of 
the linear evolution equations, which are now instead replaced by a hierarchy of nonlinear 
equations that go under the name of the Balitsky-JIMWLK equations [22-27]. With certain 
simplifying assumptions one obtains the compact Balitsky-Kovchegov (BK) equation [28] 
which can more easily be used for phenomenology. While the CGC formalism has so far 
only been written down to leading logarithmic order in x, in recent years Balitsky and 
Chirilli have derived the BK equation at the next-to-leading order accuracy as well [29,30] 
(the quark contribution was calculated earlier in [31,32]). It is therefore hoped that one 
can thus use these results to do phenomenology taking into account both the large NLL 
corrections to the evolution, and also the nonlinear corrections which are expected to be 
important at small-x. Unfortunately, however, the full NLL BK equation is extremely 
complicated and it has so far not been possible to solve it even numerically. 

We shall not present here a numerical solution of the full NLL BK equation, but rather 
use a very simple, but powerful, method to effectively take into account the nonlinear 
corrections in the full NLL BFKL evolution. We will namely solve the full NLL evolution 
using a so-called saturation boundary which allows us to extract the universal properties of 
the full nonlinear solution such as the energy dependence of the saturation scale, Qs- The 
boundary method was originally used in analytic studies of the saturation scale Qs [33,34], 
and it has since been understood that the success and justification of the method can be 
attributed, at least for a fixed coupling and at asymptotically high energies, to the formal 
correspondence between the small-x physics and the class of phenomena referred to as 
reaction-diffusion processes in statistical physics [35] (the evolution equations are said to 
exhibit "traveling wave" solutions). In [36,37], however, it was explicitly demonstrated 
via numerical solutions that the method works also for a running coupling and for non- 
asymptotic energies, and see also [38] for a recent study using the analytic methods of the 
traveling wave solutions to investigate the universal asymptotic properties of the small-x 



evolution beyond the leading order. 

When linearized, the NLL BK equation reduces exactly to the NLL BFKL equation, 
and the kernel therefore has the same eigenvalue as the NLL BFKL one [39,40]. The 
NLL BK equation thus contains the exact same double and triple poles which lead to 
instabilities in the NLL BFKL evolution. It is therefore reasonable to expect that the NLL 
BK equation will suffer from the same problems which plague the NLL BFKL evolution. 
We will here demonstrate that the NLL BFKL evolution is indeed unstable also in the 
presence of saturation effects, and thus we can conclude that the nonlinear corrections 
associated with the physics of saturation, contrary to some earlier hopes, do not cure the 
unstable linear NLL evolution. Moreover, we show here that the NLL corrections not 
associated with the running of the coupling (namely those that stem from the nonsingular 
parts of the DGLAP splitting function and the energy scale terms) are extremely important 
and that they cannot be neglected in any approximation. 

We indeed find a very strong reduction of the saturation scale, Qs{x), when the full 
NLL corrections are included. The fixed coupling results lead to the plots in figure ^ 
where we find that Q^ is reduced by around two orders of magnitude at y = In 1/x = 12 
as compared to the leading order result. Actually, the fixed coupling NLL evolution is 
highly unstable so that one can very well expect the full nonlinear evolution to be even 
more unstable. This makes it rather hard to sensibly identify the saturation scale, at least 
using the widely accepted definitions found in the literature. In this case the saturation 
boundary we apply has a stabilizing effect because the solution in the unstable momentum 
region is fixed to a certain value imposed by hand. We would certainly not expect the 
full nonlinear evolution to manifest such a regularity. Even in this case, however, we do 
find that the instability eventually kicks in, as Y grows larger, and the solution starts 
to exhibit non-sensible features. We therefore can extract the saturation scale only for a 
limited interval in Y . It should be mentioned though that we in this case have chosen the 
value Ois = 0.2 where ots = agNc/ir, with a^ the QCD coupling. It is perfectly possible 
that perhaps for a smaller, but totally unrealistic, value of a^ we could always define a 
saturation momentum; recall that for a^ < 0.05 the Pomeron intercept is positive and real, 
leading to an exponential growth of the solution. 

Part of the NLL corrections are included in the running of the coupling and when we 
allow the coupling to run, we find a somewhat milder suppression which is however still very 
strong, around a factor 7 for Y = 12. We find that the full NLL evolution is extremely 
sensitive on the precise choice of the running of the coupling. While we see that some 
choices give stable and reasonable results, other choices give very unstable results leading 
to wildly oscillating solutions for the transverse momentum distribution obtained from the 
gluon Green's function. Keeping in mind that the differences in the precise choices of the 
scale of the running coupling for the NLL kernel are formally of N^LL and N^LL order, we 
see that the evolution is extremely sensitive to the higher order corrections. Moreover it is 
also very sensitive to the minimum A; = |fc| used in the numerical implementation; lowering 
this minimum beyond some limit causes the solution to exhibit very strong oscillations 
which make it extremely unstable. 

It therefore seems that some type of resummation procedure as done in the linear 



case [11-20] is again necessary for stabilizing the evolution. Such a resummation in the 
nonlinear case is likely to be a very complicated task which is beyond the scope of the 
present paper. We will not present a full solution of the problem here, but instead we 
take the much simpler approach of studying the RG improved evolution in the presence 
of the saturation boundary. This was already studied in a semi-analytic way in [34], but 
only in the asymptotically high energy regime, and therefore potentially important pre- 
asymptotic corrections might have been missed. The application of the boundary requires 
some care with the choice of the momentum scales relevant for the process under study since 
the NLL kernel depends on the scale choice. Consequently, the resummation procedure 
also depends on the exact scale choice. The saturation boundary explicitly introduces 
an asymmetry since it acts as a cut-off on the lower values of the transverse momentum, 
modifying the linear solution asymmetrically. Similarly, the BK equation describes an 
asymmetric situation where a rather small probe, such as a virtual photon characterized by 
its virtuality, scatters off a much larger target characterized by a much smaller momentum 
scale. In this paper we consider only the asymmetric situation where a probe with a large 
scale Qa scatters off a target with lower scale Qb as is the case in Deep Inelastic Scattering 
(DIS). 

The large x terms present in the RG improved equation have a significant effect on the 
rapidity evolution. For the smallest rapidities the evolution is significantly slowed down, 
and is even negative in some k region. This behavior is due to the interplay between 
leading and non-leading terms which contribute with opposite signs. This behavior also 
manifests itself in the gluon splitting function extracted from the evolution which exhibits a 
characteristic "dip" when the splitting function is plotted as a function of the longitudinal 
momentum fraction [41]. For small and fixed coupling one can do an analytic estimate 
for the location of this dip, and one finds that it occurs when Y ~ l/y/al. As a result 
of collinear resummations, this value is, not surprisingly, parametrically far even from the 
regime where BFKL growth starts to occur, that is from Y ~ l/a^. Thus, the dip should 
have no consequences for the analytic solution when we enter the BFKL regime and in 
particular in the asymptotic one, but it is very important for phenomenology since the 
rapidity window covered by the dip region is non-negligible for realistic values of the total 
rapidity separation. We find for example that the resulting saturation scale Qs plotted 



in figure |11|, which represents the main results in this paper, stays constant, fixed by the 
initial condition, a few units in rapidity Y . Also, if we plot the resulting distribution in 
transverse momentum from the RG evolution, we find that saturation plays a smaller role 
in the evolution at these rapidities, as is manifest from the results in figure |l^ where it is 
seen that the front of the solution essentially progresses with the same speed in both the 
linear and nonlinear cases. 

The behavior of the RG improved solution can be compared with that of the pure 
LL and NLL solutions where saturation plays a bigger role. In the LL evolution, supple- 
mented with a running coupling, the nonlinear corrections are rather important and they 
significantly reduce the front velocity {i.e the rate of change with rapidity y of a point of 
fixed value for the k distribution). In the NLL case the reduction is much smaller but still 
visible even for phenomenologically relevant values of y < 15, as manifest in figure 0. As 



mentioned above the difference in tlie RG improved case is on the other hand smaller for 
the same values of Y. What this implies for the saturation corrections in the RG improved 
case is that, as already mentioned, they set in with a delay in rapidity. Needless to say, 
this behavior has interesting consequences for the search for saturation effects in experi- 
mental data where the rapidity available is rather limited. To make clear statements on 
the observed phenomenology, however, we would need to do a more careful analysis where 
the undetermined parameters and inputs in our approach are set by fitting data. 

The paper is organized as follows. In the next section we go through the BFKL formal- 
ism at both leading and next-to-leading order. We describe the choice of the asymmetric 
scale for the next-to-leading order kernel, and we outline the numerical procedure and the 
extraction of the saturation scale from the numerical solution. Then in section ^ we present 
the results of our numerical solution for both the fixed and running coupling evolutions for 
the LO and NLO evolutions in the presence of the saturation boundary. Having demon- 
strated the instability of the NLO evolution we then go on to discuss the resummation 
procedure used in our analysis in section |^. We present the exact resummed evolution 
equation which we solve, and we then present the solutions for the saturation scale and 
the Green's function, demonstrating the suppression of the saturation momentum at small 
values of the rapidity. Finally in section ^ we briefly summarize the main findings of our 
paper. 

2. NLL BFKL with the boundary 

2.1 General formulation 

Let us start this section by recalling the general formulation in QCD of the Regge limit of 
high energy scattering. Studies of 7*7* scattering lead to the formula for the total cross 
section [2] (see figure ||) 

a^'^'is, Qa, Qb) = I ^.(^y I^^ ^a{Qa, fci) G{co; k^M) ^b{Qb, ^2) , 

(2.1) 
where the functions ^A,B{Qi,kj) are the impact factors for the photons A and B with 
virtualities Qa and Qb respectively. The exact choice of the scale sq in the Mellin integral 
is arbitrary at leading logarithmic order but is important for the next-to-leading order 
calculation. The function G{uj;ki,k2) is referred to as the "BFKL Green's function" (or 
"gluon Green's function"), and should be thought of as the gauge invariant generalization 
of the vacuum expectation value of four off-shell gluons. It satisfies the BFKL equation, 
which can be written (in the case of forward scattering) [1-3] 

ujG{uj;k,ko) = 6^{k-ko) + / -^ K{k,k')G{io;k',ko) , (2.2) 

where the kernel of the equation is known up the next-to-leading logarithmic (NLL) order 
inlnl/x [4-9,42-44] 

K{ki,k2) = Ko{ki,k2) + KiikiM) +0{a\ii^)). (2.3) 



<l'j(Q_i,ti) 




Figure 1: Regge-type factorization formula for tlie cross section. The slashed gluon lines indicate 
reggeized gluons. 



It is here understood that Kq and Ki are of order q^ and a^ respectively (as clear from 
equations ( |2.7| ) and ( ^ ) below). 

We shall be solving the BFKL equation using both kernels Kq and Ki. In addition we 
will be interested in studying the nonlinear evolution equation obtained after introducing 
a saturation boundary which modifies the action of the linear kernel K. We will explain 
this procedure further below. Let us mention that we will generally be solving the BFKL 
equation with a generic initial condition in k. This corresponds to defining a new function 
T{uj, k) by 

/i2 r, 

which then satisfies the equation 

ujT{uj,k;QB) = 7^ h / — ^ K[k,k)F(uj,k ]Qb) ■ 



(2.4) 



k"^ J IT- 

This implies that we may write the cross section (|2.1D as 



(2.5) 



cT^^^is, Qa, Qb) = I^^ (^) " / ^'^a{Qa, k,)H^, fci; Qb) . (2.6) 

In the following we shall keep the dependence of J-" on Qb implicit. 

Let us now turn to the explicit expression for the BFKL kernel up to the next-to- 
leading order. We consider the solutions which are averaged over the angle. The leading 
logarithmic (LL) order kernel (after the angular averaging) is given by 



^i^o(^i,fc2)/(A;i) =a,(^2)|d^2_l_j^^(^2) _2^ 



.^11^2 ) 



f{kl)), (2.7) 



where we remind once again that as = asNc/TT. In expressing the exphcit form of the action 
of the kernel, we have above introduced an auxihary function / to simphfy the notation. 
We also use // to generically denote the scale of the strong coupling. Note that the leading 
order kernel does not have any dependence on fi since any difference in the choice of scale 
is formally of next-to-leading order. The leading order kernel is therefore scale-invariant. 
Physically this can be interpreted as the fact that in the leading order one is taking the 
limit of infinite collision energy which implies that any other finite scale in the process (for 
example transverse momenta, masses, etc) can be neglected. The next-to-leading order 
kernel is instead given by 



d^2 

TT 



^Ki{kiM)f{kl) 



, r,. 2min(A;?, fco] , r, 




Here Uf is the number of quark flavors, and Li2 is the dilogarithm function. The scale 
dependence on /i in the expression for Ki is related to the running of the QCD coupling. 
Starting from the NLL order the kernel K{ki,k2; fJ-) is thus no longer scale-invariant. 
The above form of the kernel was obtained for the so-called symmetric scale choice, see 
for example [4]. This means that the solution G{uj;k,kQ) to the BFKL equation with 



the NLL kernel above should be used for the computation of the cross section in (2.1) 
with So = QaQb- a physical example of this case is given in 7*7* scattering where the 
virtualities of the photons are comparable. On the other hand for DIS, still using the same 



formula for the cross section (|2.lD , the scale Qa is the virtuality of the photon and Qb is 
the scale characterizing the hadron target. Therefore for this situation Qa ^ Qb^ and the 
choice of scales will be asymmetric, i.e. sq = Q\- 

2.2 Scale choice in the presence of the saturation boundary 

As will be demonstrated in the numerical analysis, the scale choice in the presence of the 
saturation boundary leads to sizeable variations in the form of the solution. This is so 
because the scale choice alters the form of the kernel (2.8). In the Mellin space, it is linked 
to the fact that the scale choice changes the terms which contain triple collinear poles [4,11]. 
When combined with the nonlinear evolution this can lead to sizable differences due to the 
fact that the boundary is also asymmetric with respect to the infrared and ultraviolet 
regions (when impact parameter is not taken into account). To be precise, in the case of 



the translationally invariant BK equation the effect of the nonhnear term is such that it 
cuts off the infrared region of momenta. On the other hand, in writing down the exphcit 
form of the NLL kernel (^]^) , it is impUcitly assumed that the evolution is symmetric with 
respect to the two momentum scales, fci -f-T- A;2, which is not the case of the BK evolution. 
Therefore the correct treatment of the energy scale choice is linked with the problem of 
the symmetry of the evolution with respect to the target and projectile. This problem is 
rather difficult and it is plausible that for the complete solution one needs to take into 
account the other contributions, like Pomeron loops, which will guarantee the symmetry of 
the evolution, [45-47]. We are not going to address this important, and difficult, issue here 
but rather pick a scale choice that is relevant for the DIS process off a nucleus for which 
the BK evolution is supposed to be the correct treatment. Therefore we adopt the choice, 
So = Q'^, which is relevant for the DIS process where Q ^ Qo, with being Q the virtuality 
of the photon whereas Qq characterizes the target nucleus or a proton. Consequently one 
needs to perform the corresponding scale change in the BFKL kernel above, which amounts 
to shifting the characteristic function at the NLL level [9], 



^h) = Xiil) -2xo(7)Xo(7) • 
The functions in Mellin space are here defined as 



(2.9) 



OisXoil) 
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(2.10) 



and 



alxiil) 



dk 



^K.(fe,fe)(^^ 



2\ 7-1 



The explicit expressions for the characteristic function are 



(2.11) 



Xo(7) = 2^(l)-V(7)-V'(l-7), ^=f^ 



(2.12) 



at the level of the LL approximation, and 



XI (7) = -7j [Xo(7) + Xo(7)] - 7X0(7) - 7 
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+ U^ - T^ ^0^ + 9^3 + J-- 
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4 \^sin7r7y 3(1 — 27) 
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7(1-7) 



(l + 27)(3-27) 



E(-ir 



n=0 



V'(n + 1 + 7) - V(J0 V'(n + 2 - 7) - V(l) 



(n + 7)2 



(n + 1 -7)^ 



(2.13) 



at the NLL level. In the above equation we have set n/ = and this will be our assumption 
for the rest of this paper. Obviously, for phenomenological applications the quarks should 
be included in the evolution. For the purpose of this work, however, we will assume that 
the dynamics is constrained to the gluon sector only. 



The corresponding momentum space representation of the shift part of the kernel is 
expressed as 



/shift = 2 / r^—f{u) 
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(2.14) 



(2.15) 



(2.16) 



It can be readily checked that the above equation with 

fin) = u^-' , 

gives 

-^shift = -2xo(7)Xo(7) ■ 

2.3 Numerical implementation 

For the numerical implementation, the explicit form of the kernel ( |2.8D is not very suitable 
due to the large cancellations between the terms which constitute the contribution to the 
nonsingular part of the DGLAP splitting function. In particular there are large superlead- 
ing logarithms which cancel between terms in third and fourth lines of ( [2.81 ). To simplify 
the numerical procedure, and to obtain an accurate solution, we can instead rewrite these 
terms in a suitable way. We start by using the following form of the integral 



dx 






In 



1 



1 



1 



k^h 



1K2 



k k^ 

ln-^tan~^-^ + 29Li2 



(2.17) 



where A;<(>) = 
following form 



min(max)(A;i, A;2) to express the difficult parts of the kernel (|2.8D in the 
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rf\ A^^ t\i. 
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22 ]_f^ 4 
^+32"32V^ + ^ 



(2.18) 
''> "'</!['*'< V '*'> / V '^> J . 

The functions Si and S2 are series expansions of (tan~"'^(x) — x) and (9Li2(x) — x) re- 
spectively. We have checked that it is sufficient to retain only around twenty terms in the 
expansions to get accurate results. 

2.4 Saturation scale from the boundary method 

One of the main objectives in this paper is to extract the saturation scale Qg from the NLL 
and RG improved evolutions using the boundary method. In this section we therefore first 
describe the definition of the saturation scale used, and having done that we then describe 
the precise numerical method in which the definition is employed to obtain Qg. 



2.4.1 The definition of the saturation scale 

The usual definition of the saturation scale follows from the solution to the nonlinear 
BK equation. The object satisfying the BK equation is the coordinate space scattering 
amplitude^ M{s, r) for a dipole of size r. Here we denote s/sq by s. In this case one can 
define the saturation momentum Qs as the scale which separates the regions where N is 
"small" and thus follows a linear evolution, and where it is nearly saturated at A/" = 1 
and its evolution is completely nonlinear. The exact definition of Qs is always somewhat 
ambiguous since it depends on the precise value of N where one chooses to separate the 
linear and the nonlinear regions. 

In determining Qs one defines first the critical dipole size rs{s) by N{s,rs{s)) = c < 
1 where c is a given constant smaller than 1. Then one can take Qs = C r~^ where 
C is another constant. In the Golec-Biernat-Wusthoff model [48, 49] for instance, one 
chooses C = 2. The constant c can be chosen to be around 0.1 — 0.5, the exact value 
will determine only the normalization of Qs which is never under full control theoretically. 
What is determined by the perturbative evolution is the scale dependence of Qs, i.e. the s 
dependence. Note that the saturation of AA does not necessarily imply the saturation of the 
cross section (t(s, r) of the dipole. The saturation of the cross section is a nonperturbative 
problem which is not solved by the perturbative nonlinear evolution equation for M. To 
obtain the cross section a one needs to integrate M over the impact parameter b 

aA^{s,r) = 2 j d%N{s,r,b). (2.19) 

Obviously the behavior of N at large h is nonperturbative and needs to be modeled phe- 
nomenologically. In the GBW model the dipole cross section a is taken to be 

adip(x,r) = aoM-exp(-^-^n . (2.20) 

Here one sets s = x"^. In this case r^ is the length scale above which the cross section of 
the dipole saturates to the constant ctq. If one assumes that the h dependence factorizes 
as J\f{r,h) = J\f{r)S{h) then uq = 2Jd%S{h), and in (|3^) clearly r = rs{x) is a line of 
constant N . Of course the assumption that ctq does not evolve with x implies in this case 
the complete saturation of the total cross section. It is well known that the x dependence of 
the fitted rs{x) agrees in form with the one obtained from the solution to the BK equation, 
using the definition of constant N , but also that the leading order evolution gives a much 
too steep growth with \/x. 

In this paper we are, however, solving the momentum space BFKL equation using the 
saturation boundary. The question then is how exactly we should define the saturation 
scale from the solution to the evolution equation. It is possible to consider different choices. 
For example, in [50] Qs was defined as a line along which the difference between the linear 
and the nonlinear solutions to T (in that case obtained from the leading order BK equation) 



^Equivalently the BK equation can also be written for the dipole "S-matrix" defined as S = 1 — A/", see 



equations (3.8) and (3.9) 



10 



is of a certain magnitude. This, however, gives an energy dependence somewhat different 
than the coordinate space definition. 

We shall therefore choose a simpler prescription whereby just as in the coordinate 
space (similarly to (|2.2[l|) ) we define Qs via 

J"(s,(5s(s)) = const X ctq • (2-2f) 

Notice that according to ( |2.6| ), F equals the cross section when the impact factor <1>^/A;^ is a 
delta function in the transverse momentum. That would be the case if the scattering object 
A itself is a parton. Thus one can think of J- as the cross section of a parton impinging 



on the target particle. Therefore the definition (2.21) bears a certain resemblance to the 



definition ( 2.20] ) . Another partial motivation for this simple choice comes from the fact that 



this definition gives an x dependence of Qs which is rather consistent with the empirical 
dependence extracted from data. We are here not concerned with the exact value of o"o 
which we shall leave unspecified^. The exact value of the right hand side of ( p.21 ) will 
determine only the absolute normalization of Qs- 



Another motivation that can be given for the definition ( |2.21 ) directly in momentum 



space is, at least quasi-classically, a number density in the transverse phase space. In the 
Color Glass Condensate model, the equivalent number density 

"" (2^22) 



of the classical fields saturates when it is of order l/us (see for example [51]), and it leads 



to a definition of the saturation scale exactly as in (2.21) when integrated over the impact 
parameter. Now, our T, which satisfies the BFKL equation and which can be used in a 
formula like ( p. 61) to calculate the total cross section, cannot literally be thought of as a 
number density in the phase space. It rather corresponds to a cross section at the parton 
level which is conceptually a different object than a phase space number density, since the 
gluon Green's function does not have an operator definition which would match exactly 
the definition of a number density. 

The object defined as a phase space number density, which is essentially the expectation 
value of the field-strength tensor (F+*F+*) in light-cone gauge (see equation (2.15) in [51]) 
does, however, also satisfy the BFKL equation at the linear level (see (3.59) in [51]). On the 
other hand it does not satisfy the BK equation at the nonlinear level, but some equation 
more complicated than BK (this equation can be obtained by applying the JIMWLK 
kernel to the operator definition of the quasi-classical number density given by equation 
(2.18) in [51]). The point now is that it does not matter for our analysis what the exact 
nonlinear equation is, no matter how complicated it may be. The boundary method is 
a generic method which can extract the universal properties of the full solution, such as 
the s dependence of Qs, and it therefore works for nonlinear equations which are in the 
universality class of the BK equation, that is equations whose linear parts are driven by 



^As we are here not doing any phenomenological fit this wiU not be important. The exact value cto 
would be determined as a fit parameter in any practical application. 
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the BFKL kernel. Therefore if we simply regard J- only by its property that it satisfies 
the BFKL equation, that is to say if we simply forget the would be operator definition of 
the gluon Green's function and only consider the equation it obeys, then we do not have 
to care whether it actually corresponds to the phase space number density or to the dipole 
scattering cross section, we will obtain the same solution for the saturation scale. Since 



the definition (2.21) works in both cases we conclude that it is indeed the optimal choice 
for our problem. 

We therefore now move on to the exact implementation of the boundary which effec- 
tively takes into account the missing nonlinear effects. 

2.4.2 Application of the saturation boundary 

In accordance with the definition (2.21[), we shall apply the saturation boundary in our 



numerical treatment as follows. First of all, as explained in the previous section we need 
not specify the exact value of (Tq for the purposes of our study. Let us therefore define the 
new function F(s, k) by 

F(|,fc) = :^(!l^, (2.23) 



which means that the condition ( 2.21 ) defining the saturation scale now simply reads 



F(s, k = Qs) = const. (2.24) 

Obviously F satisfies the exact same equation as J-. 

We then define first a so-called "critical value" , c, which for convenience will be taken 
as a number of order 1. This number can be close to the constant in ( |2.21] ) but ideally it 
should be slightly smaller^. At each step in the numerical solution of F we then define the 
corresponding critical transverse momentum scale kc via 

F{s,kc{s))=c. (2.25) 

The boundary will now be applied to those transverse momenta which are below the critical 
scale by some magnitude determined by a second parameter A. More precisely, if by p 
we denote the logarithmic scale p = ln(/c^/A;Q) then F{s,p) is forced to satisfy the given 
boundary condition for all 

p<Pc-A, pc = \xi{kl{s)/kl) . (2.26) 

We shall choose the arbitrary /cq as the minimum k of our numerical computation. 

In the original analysis of the method in [33] the boundary was chosen to be totally 
absorptive, i.e. 

F(s, /)) = for all p < p^s) - A . (2.27) 

This is the choice most appropriate for the analytic analysis since, after the leading expo- 
nential behavior has been factored out, the problem can then be formulated as a random 



^This is so since our c is loosely speaking the value where the effects of saturation first starts to play a 
role while Qs can be thought of the scale below which the evolution is really dominated by the nonlinearities. 
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walk in the presence of an absorptive wall [33]. For the numerical solution on the other 
hand we can chose any condition which cuts off the power-like growth of the linear evolu- 
tion with s. Actually the totally absorptive boundary has to be applied with some care, 
in the numerical treatment, when the evolution is not "fast enough". What happens in 
this case is that when all the contribution below the boundary Pc{s) — A is cut completely 
by hand, the next step of the evolution may not be strong enough to push the solution 
above the critical value c. If this happens, then at this next step the boundary does not 
get applied by definition, since the solution completely falls below c. As saturation then 
suddenly switches off, there is an accumulation of the solution just around the boundary 
which can give a "spike" in the solution (that is the solution suddenly jumps high above 
the critical value before it is set to zero), and moreover some of the previous points in k 
where the solution was set to zero may now grow which in turn implies that the front may 
actually, in a single step of the evolution, move in the wrong direction (to smaller k). We 
have observed this somewhat peculiar behavior also in the case of the CCFM evolution 
equation where the evolution is suppressed at large k due to the restriction imposed by 
angular ordering. It also appears to be the case in the RG improved evolution we study 
below. On the other the leading order BFKL evolution always seems fast enough so that 
this problem never appears. We do not find that this is a major obstacle for the numerical 
implementation of the totally absorptive boundary, but it does require some care in the 
precise treatment. We therefore also consider a second boundary which is more straight- 
forward to implement numerically. In this case we simply freeze F at the boundary, and 
one can show that the two implementations are equivalent in the asymptotic regime. That 
is, at each step in the solution we let 

F{s, p) = F{s, pc - A) for ah p < pc{s) - A . (2.28) 

Let us here emphasize that these choices should not be interpreted literally as to how 
exactly saturation would act on F. For example, freezing F at a constant value does 
not imply that one should think of the dynamics as representing the physical result of 
the saturation of some gluon occupation number at a fixed value. The whole idea of the 
boundary method is that it does not matter how the linear growth is cut off. We have 
chosen these two boundaries for their simplicity, not because they would represent a more 
accurate representation of the true nonlinear terms compared to other possible choices. 

Studying the evolution in the presence of these two boundary conditions we shall ex- 
tract the universal properties of the solution that are independent of the precise conditions. 
As is evident from our construction, c and A can be thought of as free parameters related 
to the freedom in choosing the precise way via which the missing nonlinear terms are ac- 
counted for. In reality, however, these parameters are not completely free. First, they 
are correlated as A ~ Inl/c [33]. Secondly, the value of c in (2.25|) cannot be completely 



arbitrary. As we noted above, c should be of the same order of the constant used in ( 2.21 ), 
and the latter is of order 1. We have chosen our default value to be c = 0.4, but as in the 
previous works [36, 37, 52] we have also considered the possible sensitivity of the solution 
to different choices of c (and A). 
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Figure 2: Solutions to tlie leading and next-to-leading BFKL equations for a fixed coupling 
ds = 0.2 and for Y = 2,6, 10, 14: NLL BFKL with asymmetric scale choice (sohd red), NLL BFKL 
with symmetric scale choice (solid blue), LL BFKL (dotted blue). 



3. Results for NLL BFKL with and without saturation 



3.1 Results with fixed coupling 

In this section we present the results obtained using a fixed coupling, a^ = 0.2. We start by 
studying the linear evolution equations. For the results shown in figure |2|, and in subsequent 
plots, we have chosen the initial condition 



Fo^k) = c ■ exp 






(3.1) 



with kin = 1 GeV. Here c is the parameter in ( 2.25 ) which determines the critical value 
beyond which the boundary is applied to F. We also use the logarithmic variable Y = 
In s = In s/Q"^ to denote the energy dependence of the solution. The comparison between 
the LL and the NLL BFKL solutions is shown in figure ^. Here the NLL solutions are 
shown for both the symmetric and the asymmetric scale choices. As well known, see for 
example [9,14,53,54], the NLL evolution is significantly slower than the leading order 
evolution. We also clearly see the effects of the negative Mellin space poles in the collinear 
and anti-collinear limits. For the symmetric solution (the solid blue curves in figure ^ 
the collinear pole is clearly visible in the plot as the solutions turn rapidly negative at 
moderately high values of the momentum k. For the asymmetric scale choice the shift in 
the characteristic function (|2.9| ) removes the collinear triple pole while there is still the 
double pole. As a consequence it turns negative "later" (i.e. at higher values of A;) than 
the symmetric solution. On the other hand the pole at the anti-collinear end causes the 
solution to turn negative at the smaller k values as clearly visible in the figure. 

The apparent instability of the solution suggests that the full-linear solution might 
very well be even more unstable. The precise behavior will of course depend on the exact 
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shape of the nonhnear terms and how these act, but it is clear that the existence of a 
nonhnear damping term F'^ or F"^ can potentially cause further instabilities when F turns 
negative. It is therefore not all clear whether one would obtain any sensible results. Of 
course one can imagine that the results are stable for a smaller coupling. For example it is 
known that for as ^ 0.05 the Pomeron intercept is positive and real in the NLL case, and 
this leads to an exponential growth of the solution. Such small values of a^ are of course 
hardly realistic. 

The boundary method on the other hand is very stable by construction since the 
solution below the boundary is set to a fixed value by hand. Despite this, however, we 
will see below that the final result is nevertheless unstable. Note also that, the successful 
implementation of the boundary method should be independent of the precise value of the 
critical value c. For the LL solution in figure ^ it is clear that no matter what c is chosen, 
the linear solution eventually reaches this point, and the saturation boundary therefore gets 
implemented. The difference between various values of c is simply in the normalization of 
Qs which of course is a measure of the strength of saturation (and therefore of the value of 
c). For the asymmetric NLL solution (solid red lines) in figure |2|, it is, however, clear that 
the peak of the solution is bounded due to the solution turning negative at smaller k. In 
this case if we would choose c to be say larger than 1, then it appears that the boundary 
would never get applied, and consequently it would have no effect whatsoever on the linear 
solution. Clearly this is a peculiarity of the unstable NLL evolution, and the results for 
the fixed coupling case should therefore not be taken too seriously. What we can learn 
from the fixed coupling case is on the other hand that the full non linear solution might 
very well not give any sensible results (it would of course be highly desirable if this can be 
checked explicitly). 

Keeping these points in mind we now apply the saturation boundary to the linear 
solutions shown in figure ^. In accordance with the discussion above we study the nonlinear 
evolution at the next-to-leading order only with the asymmetric scale choice obtained by 
applying the shifts in ( |2.9| ) , ( |2.10| ) and ( p.llD . In figure ^ we compare the fixed coupling 
NLL solutions with and without the absorptive saturation boundary. We see that the 
instability of the NLL solution at lower k is not removed completely by the condition of 
the absorptive boundary. For the solution at y = 14 we see that F turns negative before 
it is set to by the boundary condition. As in the linear evolution the solution also turns 
negative at higher k which is of course expected since the nonlinearities do not cure the 
unstable high-fc behavior. 

Note also that the linear and nonlinear solutions differ slightly even in the very large k 
region. It seems that this is due to the fact that the nonlinear solution with the saturation 
boundary kills the contributions at small k which via the NLL kernel contribute negatively 
to the high-A; part of the solution. If this is indeed the case, we would for the second 
boundary condition (2.28) expect that the difference compared to the linear solution at 



larger k is somewhat smaller since in that case the contributions below the critical point 



are not set to zero. We show the results for the boundary ( 2.28| ) in the left plot of figure 
^ where we can see that this is indeed the case. In the right plot we instead compare the 
solutions obtained by the two boundary conditions. Note that when Y is large enough. 
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Figure 3: Comparison of the NLL BFKL solution with (solid red) and without (dotted blue) the 



absorptive saturation boundary (|2.27D for fixed as = 0.2, and F = 2, 6, 10, 14. 
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Figure 4: Left: Comparison of the NLL BFKL solution with (solid red) and without (dotted 
blue) the frozen saturation boundary ( 2.28 ) for Y — 2,6,10,14. Right: Comparison of the NLL 
solution with the frozen (solid red) saturation boundary ( 2.28| ) versus the absorptive (dotted blue) 
saturation boundary for F = 2, 6, 10, 14, 20. Simulations done for the fixed coupling ctg = 0.2. 



also the solution obtained using the boundary (|2.28 ) turns negative at smaller k. Thus we 
see that the nonlinearities associated with the mechanism of gluon saturation do not cure 
the unstable low-/c behavior completely, even when they by construction stabilize the low 
k region. 

3.1.1 The NLL saturation scale for a fixed coupling 

As is clear from above it is hard to define the saturation scale Qs beyond a certain value 
of Y because of the severe instabilities of the solution. The results for the saturation 
scale presented for the fixed coupling case should therefore not be taken too seriously, but 
we here want to demonstrate the very large effects of the higher order corrections on the 
solution. We again mention that the full nonlinear equation might be even more unstable 
so it is not clear whether the standard notion of the saturation momentum even makes 
sense for the chosen value of 0.3- 



16 



lU 






1 




1 




1 




1 




1 


• 


1 : 
























• 




" 






















• 






■ 


^10' 


- 
















9 








- 






























eg 




























> 
















• 










- 


0) 


























- 


O 














• 












- 






























-^"10^ 


- 










• 














~i 


o 








• 


• 






• 


• 


• 


• 


• 


• i 


10° 


_ 


• 


• 


# 


• 


• 


• 












_ 


<i 


• 


• 
1 




1 




1 




1 




1 




1 



10 12 14 



Figure 5: The saturation scale Q^{Y) extracted from the evolution using the boundary ( 2.28 ), for 
LL BFKL (red, upper circles) and NLL BFKL (blue, lower circles), with fixed coupling a^ — 0.2. 



The results are shown in figure |5| where we have extracted the saturation scale, Qs{Y), 
from the leading and next-to-leading order solutions using the boundary ( p. 281 ). While 
the leading order evolution gives a very rapid, exponential, increase of QsiY) with Y, 
the next-to- leading order evolution leads to a strongly suppressed result. Thus we can 
expect the NLL corrections to the BFKL kernel to have dramatic effects for the study of 
the nonlinear evolution as well, but the analysis cannot be complete until the important 
running coupling effects are taken into account. As we will see below, the running of the 
coupling has a rather large effect on the evolution and on the result presented above. 

3.2 Results with running coupling 

In the case of the running coupling we are faced with the question as to how exactly choose 
the scale of the coupling. For the scale dependence of a^ we shall use the one loop result 
resummed to all orders, that is 



as{k) 



1 



Mn((A;2 + ^2^)/A2) 



(3.2) 



where we have inserted an infrared regulator of the Landau pole. We shall by default set 
fJ'iR = 0.7 GeV. We should immediately note that the linear evolution is sensitive to this 
parameter since a smaller /xir implies an enhanced contribution from smaller momenta 
which generally speeds up the growth. To check the sensitivity to ^m we have also run the 
simulations with //ir = 0.4 GeV, which, as expected, speeds up the growth of the linear 
solutions, but we find that to a very good accuracy it does not affect the nonlinear solutions 
obtained from the saturation boundary which are therefore rather robust with respect to 
the regulator. 

As can be seen from the NLL BFKL equation, the natural scale in the leading part 
of the kernel is given by the transverse momentum of the real gluon, q. The choice in the 
NLL part of the kernel is on the other hand rather arbitrary since any difference in the 
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scale choice is formally of N^LL and N'^LL order. However, as we shall see, this formally 
higher order difference is extremely large and the solution therefore very sensitive to the 
exact choice. 

We have investigated different prescriptions for the running coupling. Let us write the 
BFKL kernel as in ( p.3| ) but this time extracting out the factors of as- 

K{k, k') = as Ko{k, k') + al Ki{k, k'). (3.3) 

In studying the running coupling NLL evolution we have considered the following choices 
(/;;> = max(A;, k') as before) 

A: 



B 
C 
D 



as{q^) Ko{k, k') + al{kl) Ki{k, k') (3.4) 

as{kl) Ko(k, k') + alikl) Ki{k, k') (3.5) 

Q,(A:2) Koik, k') + al{k^) Ki{k, k') (3.6) 

as{q^) Ko{k, k') + al{q^) Ki{k, k') (3.7) 

where in choice D, the real momentum q is used in all the real terms only (the virtual terms 
are diagonal in k). In all cases the kernel Ki has been adjusted so that the expressions are 
the same at the NLL level. 

We find that the choices A and B give rather well behaved solutions but that C and 
D lead to very unstable results which rapidly turn negative and oscillate over very large k 
intervals. The differences between the choices A and B are not that large and so here we 
will only present results obtained from choice A. The fact that this choice gives a stable 
result is consistent with the findings of [14] where the different scale choices of the running 
coupling were also investigated. We note that in this choice the NLL kernel does not 
contain any terms which depend on the beta function coefficient. 

The results extracted using the scale choice A above are shown in figures |^ both for the 
linear and the nonlinear solutions. For the linear case we again find a very large difference 
between the leading and the next-to- leading order solutions. The leading order solutions 
are here obtained using a running coupling which runs with the scale q as clear from choice 
A in ( |3.4| ). In this case the NLL solution seems better behaved than in the fixed coupling 
case studied in the previous section. It should, however, be kept in mind that while choice 
A is stable, other choices like C and D lead to unstable results. Thus the natural question 
arises whether the NLL evolution including the running coupling has any predictive power 
since it is extremely sensitive to corrections which are formally of higher order. 

In the right plot in figure |^ we show the nonlinear solutions obtained after applying 
the frozen boundary condition ( |2.28| ). Apart from the solution at the lowest Y , we see that 
the running coupling leading order evolution (rcLL) has essentially the same slope as the 
full NLL solution but that it again grows rather more rapidly with Y . We also compare 
directly the linear and nonlinear solutions obtained from the boundary in figure 0. In 
this plot, the differences between the solutions with saturation and the linear solutions 
are better visible. In the leading logarithmic case, application of the boundary affects the 
region far away from the boundary much more than in the next-to-leading scenario. This 
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Figure 6: Full NLL BFKL solution (solid red) versus leading order running coupling solution 



(dotted blue) with a running coupling using the scale choice A in (3.4) for Y — 2,6, 10, 14. Solid 
black curve is the initial condition at y = 0. Left: Linear evolution. Right: Nonlinear evolution 



using the boundary (2.2g). 
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Figure 7: Linear BFKL solution (dotted blue) versus solution with the boundary (solid red) with 



a running coupling using the scale choice A in (B.4) for Y ^ 2,6, 10, 14. Left: LL case. Right: NLL 



case ( |2.28| ). 



has a prominent effect on the speeds of the front evolution which in the next-to-leading 
case is much less affected by the saturation corrections. 

The large sensitivity to the different scale choices in the running coupling is not un- 
expected. Given the large values of the running coupling the truncation of the small-x 
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perturbative expansion will lead to large uncertainties. These could be avoided by setting 
the scale choice via the BLM [55] scheme which reduces these uncertainties. Indeed, it has 
been demonstrated that the NLL small-x evolution can be stabilized in this scheme [56]. 
We note however that the resulting scale in the BLM scheme has a rather large numerical 
coefficient which reflects substantial differences between the MS scheme and the BLM- 
MOM scheme adopted in [56]. A related observation was done in [57] where the principle 
of minimal sensitivity was applied to the NLL BFKL in the process of electroproduction 
of two vector mesons. There it was found that the optimal choice of the renormalization 
scale is about ~ lOQ with Q^ being the virtuality of the colliding photons [57,58]. One 
could explain this unnatural choice of scale by the fact that the higher order subleading 
corrections are effectively taken into account via this procedure. This was later confirmed 
by redoing the analysis [59] using the collinearly improved, resummed kernel [18], in which 
case the resulting optimal scale turns out to be ~ 3(5, which is much closer to the typical 
scales involved in the process. What all this shows is that the higher order corrections are 
generally very large, and that the NLL evolution effectively stabilizes only by choosing a 
scale of the running coupling which is much larger than the natural choice dictated by the 
relevant physical scales in the process^. 

We also should add that the solutions in general are very sensitive to the lower cutoff 
on momentum k. In the simulations presented in this section we used the cutoff fe^jn = 
0.1 GeV . This sensitivity is of course due to the large value of the coupling in this regime. 
It is worth noting, however, that the NLL solution does show some instability with respect 
to the variation of this cutoff. Below /c^j^^ = 0.05 GeV we find that the low k part of the 
solution turns negative and then oscillates very strongly leading to severe instabilities (also 
for the otherwise stable choices A and B\ While also the value of the LL solution (and 
resummed which we analyze later) increases with decreasing cutoff, in that case there is 
no instability of the solution as in the NLL case. We have therefore chosen the femin here 
so that the NLL evolution gives stable results. 

3.2.1 The NLL saturation scale for a running coupling 

We easily extract the saturation scale using the definition ( p. 24 ) from the solution shown 



in figure |6[ The results are shown in figure ^ where the running coupling leading order 
results are compared with the full NLL results. As we can see the differences between the 
leading order and the next-to-leading order kernels are smaller than in the fixed coupling 
case, but we also note that they are again generally rather large. At y = 2, Qg differs by 
around a factor 1.5 between the running coupling leading result and the full NLL result, 
while at y = 18 this difference has grown to a factor of almost 10. Thus we clearly see 
that the next-to-leading order corrections beyond that of the running of the coupling are 
very important and cannot be neglected. 



Moreover, as mentioned above the scale choice A (3^) by which the results in figure 



|8| have been obtained is actually the more stable one giving a steady growth of F with Y . 

*It is indeed not strange that the NLL evolution would be well behaved in this case. Since the scale of the 
coupling is so large, and consequently its strength so small, the higher order corrections are automatically 
suppressed and presumably not important. 
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Figure 8: The saturation scale QHY) extracted using tlic boundary ( 2.2§| ) for the running 
coupling solutions to the LL kernel (red, upper circles) and the full NLL kernel (blue, lower circles). 



While the running coupling leading order solution does not give a markedly different Qs 
for other choices of the scale of the coupling, the NLL solution changes dramatically, so 
that we cannot even sensibly extract a saturation scale. 

One might have already wondered what is the origin of these pathologies in the NLL 
evolution. To this end, it is instructive to look at the the NLL nonlinear equation in 
coordinate space, that is the NLL BK equation, for reasons that will be clear in what 
follows. The BK equation at leading order reads 



dS.. 



xy 



AY 



/ d Z JVlxyz\>^xz>^zy ^xy) ; 



(3.8) 



and has a straightforward interpretation. Sxy is the S'-matrix for the scattering of a color 
dipole {x,y) of a generic hadronic target. Under an increment dY in rapidity, a soft gluon 
is emitted at the point z and we view it as a quark-antiquark pair at large Nc- Thus the 
parent projectile dipole can split into two, {x,z) and {z,y), which subsequently scatter off 
the target as suggested by the first term in the r.h.s., with the second corresponding to a 
self-energy correction. The probability for the splitting is of order 0{as) and is given by 
the kernel Mxyz which is known and is positive for any value of its arguments. 
At NLL order the BK equation becomes (for example cf. Eq.(104) in [40]) 



dS: 



xy 



dY 



d^zM^xlz 



{^xz^zy 



S: 



xy 



+ 



d^zd^wM'iizw 



{S. 



xz^zujJwy 



W 



(3.9) 



where, for the sake of clarity, we have kept only the dominant terms at large Nc- Then the 
above again has a nice interpretation: the parent projectile dipole {x,y) splits either into 
two dipoles {x,z) and {z,y) or into three {x,z), {z,w) and {w,y), which then scatter off 
the hadronic target. The probabilities for these splittings to happen are given by the two 
kernels A4^^' and Ai^^', with Ad^^' containing the leading piece A4^^' of order 0{as) plus 
the NLL contribution of order 0{a'^), while A4^^' is of order 0(a'^). 



21 



The problem is that a direct inspection of the two kernels reveals that they can both 
become negative and, moreover, large in magnitude. This will happen in the collinear 
and/or anticollinear limit, that is, when the emitted gluons at z and w are emitted very 
close to one of the parent color sources at a:; or ?/, or very far from them. Since collinear 
splittings are related to the ultraviolet behavior of the evolution, we cannot expect satu- 
ration to cure this problem of large negative probabilities. We do not analyze this in more 
detail here, since in the next section we shall see how these pathologies arise in Mellin 
space. 

Thus our main conclusion studying the NLL evolution in the presence of saturation 
effects is that the instabilities of the next-to-leading order BFKL evolution are not cured 
by the presence of the nonlinear corrections. In order to obtain sensible results out of the 
evolution equations we must therefore consider an improvement that can cure the insta- 
bilities inherent in the formalism. As a specific model we here consider the resummation 
technique presented in [14]. We thus now turn to the renormalization group improved 
BFKL evolution which we shall subsequently study in the presence of saturation effects. 

4. Resummed BFKL with the boundary 

We begin this section by describing the resummation model presented in [14]. As we have 
seen above, a resummation is needed in order to control the large higher order corrections 
in the small-x evolution. After describing the method we use, and the precise equation 
that it implies, we will go on to apply the saturation boundary and extract the resulting 
saturation scale. 

4.1 Construction of resummed kernel 

Let us first briefly recall the principles of the resummation procedure originally presented 
in [14], where more details can be found. Similar approaches to resummation have been 
developed in [15-17] with consistent results. The procedure consists of the resummation 
of the collinear singularities at NLL small-x evolution, as well as in the incorporation of 
the running coupling. The eigenvalue ( p. 13 ) contains double and triple collinear poles in 



Mellin space, which are numerically very large corrections. Its approximate behavior near 
the 7 — >• and 7 — )■ 1 poles is 

^1 ^^^ 273 2(1-7)3 72 ^ (1-7)2 ^ ^ 

Here ^i(O) = —11/12 is coming from the nonsingular part (in momentum fraction z) of 
the leading order DGLAP splitting function 

Iggi^) = as- +asAi{uj) , (4.2) 

UJ 

where Mellin transform has been defined as 

Pgg{z)z'^dz = -fggiu) . (4.3) 
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We shall also define 



^ a.. 



Pgg{z)z^dz = asAiiu) = ^gg{uj) - -^ . (4.4) 

^ 



In (4.1), the double poles come from the DGLAP splitting function as mentioned 



above, whereas the triple collinear poles come from the choice of the energy scales. The 



above collinear kernel (4.1) is negative and is responsible for 90% of the corrections of 
the NLL kernel in the case of the fixed coupling, that is when we set 6—7-0. In order to 
smooth out the behavior near the collinear poles, the resummation procedure has to be 
taken into account. The two basic ingredients of this resummation are the subtraction of 
the above collinear poles and their replacement by the resummed expression which is due 
to the nonsingular part of DGLAP, and the kinematical constraint. That is, the resummed 
kernel reads 



with 



Xresum(7,'^) = Xo(7,w) + Xcoll(7, '^) + "sXiIt) > (4-5) 



Xo(7, w) = 2^(1) - ^'(t + I) - ^(1 - 7 + f ) 

ujAiiuj) ujAi{uj) 
Xcoii(7,'^) = — -Tr + 



7 + f 1-7 + f ' (4.6) 

Xi(7) = Xi(7) + 0X0(7)— ^7— T - Xo(7): 



2 sin^(7r7) 7(1 — 7) 

where the last line contains subtractions of the triple (second term) and double poles (third 



term) which numerically coincide with (4.1), modulo the b dependent term. In contrast to 
the LL and NLL kernels, one of the important features of the resummed equation is that it 
satisfies energy conservation. Note that now the kernel eigenvalue is a function with both 
7 and uj dependence. This is because the dependence on the coupling constant for some 
terms has been traded off for the dependence on w, due to the resummation. 

The above eigenvalue is strictly correct for the fixed coupling case. In the running 
coupling case, it is better to consider the momentum representation and perform the re- 
summation of the coupling into the argument of the coupling in front of the leading order 
kernel. Based on the resummed expression in Mellin space, the following proposal was 
made for the kernel in momentum space (we use again the notation / to represent the 
action of the kernel as in ( ^.Tj) and ( |2.8D ) : 



/"'^ fdk'^]C{z;k,k')f(-,k'^ 
Jx z J \z , 

- j ^fdk'^[as{q^)K^%z;k,k') + as{kl)K^^{z-,k,k') + al{kl)Kiik,k')]f(^^ . 



(4.7) 



The different terms are as follows: 
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LO BFKL with running coupling and consistency constraint {q = k — k') 
^^ dz 



dk" 



dz f d^q 



as{q')K^%z;k,k') f{-,k' 



oisiq^ 



f(^,\k + q\\e(--k'j@ik'-kz)-Qik-q)f 



(4. 



Nonsingular DGLAP terms with consistency constraint 
"^ dz 



dk'^asikl)K^%z;k,k')f[-,k' 



X 

^z ' 
^ dz /■'=' dA;'^ .,^. k ~ ( k\ J X 



(kzy 



fc2 



k 



^s{k')z-P,Az-]f{-,k' 



z 



dz r^^/^)'dk'^ _ ,2. k'~f k'\/x ,, 



(4. 



NLL part of the BFKL with subtractions included 



d^ f "'2 ^,2„2^r> /,. U^rfx ,j\ _ 1 f^dz I ^ .2 ^2/, 2 



— / dk'^ ai{ki)Ki{k,k')f{ -,k' 



^ ,^- I dk'^ ^HK) 



67 ttA 1 
'9 ~ Y ) \k'^ - k^\ 



n\.^' 



2^< -ji^'-^k^ 



+ 



"32U'2 + A;2 + U'2 A;2J^''U/2 



A;'2 + A;2 y^ 

AU2(l-k%/k 



-4:Ai{0)sgn{k'^ -k'"]l t^ In 



3 + 



" ^VA;2 
32k'^k^ 



\k" - fc2 



+ 
1 



2 ^ 



Ifc' 



/2 



A;2| 



A:' 
dy 



■In. 



|fc'2-A:2| 



A;2 + y2fc/2 



In 



fc2 
1 + y 



i-y 



+ 



1 ^^' .T- ^4 



f -,k' 

z 



6C(3) f'dz , , 






(4.10) 



The fact that the shifts of the collinear poles are symmetric in ( [4.6| ) is reflected by the 
symmetric form of the kinematical constraint in ( [4.8D and ( [4.9D . This in turn is related to 
the symmetric scale choice, sq = QaQb- As mentioned above, in our calculation we actually 
use the asymmetric scale choice, sq = Q\. In addition, in the kernel written above one 
needs to perform additional subtractions. This is because the kernel still contains some 
residual single poles. They contribute to a residual 2-loop anomalous dimension which 
needs to be subtracted. For the purpose of this analysis we choose the scheme called 
scheme B in [14]. It consists of a modification which adds a term with the shifted pole to 
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the NLL kernel with the w-dependent coefficient 



xih)^xih) = xih) 

where 



1 1 

- + 



7 1-7 



C(0) + 



1 1 

+ 



7 + 1 1 + 



C{iv) 



Aiiu) , V(l + w)-^(l) 



w + l 



+ 



UJ 



vr 



(4.11) 



(4.12) 



C(0) = --^i(0). 



This scheme also satisfies the energy-momentum sum rule for the extracted resummed 
anomalous dimension 'jgg. 

The change in the resummed kernel in (x, fe^) space corresponding to scheme B is 
obtained by taking the inverse Mellin transform of (|4.11|), and is given by 



dz 



dk'^\^as{q^)K^^z; k, k') + as{kl)K^''{z- k, k') + al{kl)ki{k, k')^f(-, k' 



dk" 

(kz)2 k'^ 



■?.■?. k f k\ fx \ /■('=/^) dA;'^ 2/,/2n k' ( k'\ fx , 



(4,13) 



with the function S2{z) given by 



^^^")=ii;r^2+^^^'+" 



-541 + 2477^ + 72z(l + 3z) 



1441n( -1 + - j Inj- 



+ 12 ln(l - z) -1 - 2^(23 + z(-15 + 8z)) - 12(1 + z) ln(l - z) 
+ 12zln( -1 + - J Inj - J + 2z[l + z(-21 + 5z) - 61n(l - z)] \n{z) - 6(-l + 2z) \n^{z) j 



+ 144(-1 + z) 



Li,fz) + iln( i)ln 



2 \z 



{l-z) 



IT 

~6 



144(1 + 2z)Li2(l - z) ^. (4.14) 



4.2 Numerical results 



In this section we present the numerical results for the evolution using the resummed 
kernel, also including the saturation boundary implemented via the conditions ( |2.27| ) and 
(2.28) as before. Before going on to the study of the resummed kernel in the presence 
of the saturation boundary, however, let us first recall the comparison of the solutions 
for the linear LL, NLL and resummed kernels which was performed in [14]. In order to 
compare our different solutions to the gluon Green's function to those obtained in [14] , we 
use the symmetric scale and an initial distribution of the form of a discrete delta function 
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Figure 9: The gluon Green's function solutions to the hnear BFKL equation, obtained using 
a symmetric scale choice and a discrete delta initial condition. The red dashed line is the LL 
approximation, the black dashed is NLL, and the solid blue is the resummed prescription. In each 
case (including the LL case) the strong coupling constant is running as as{q) in front of the leading 
order term, and as as(max(fc, fc')) in front of the subleading terms. The scale fco = 20 GeV; all the 
solutions correspond to rapidity Y — 10. 



in accordance with ( |2.2D (rather than the initial condition (p.l|)). These results are shown 
in figure ^. In all cases as{q) is taken as the choice for the running coupling in the leading 
term and as{in.ax(k, k')) for the subleading terms (we recall that this corresponds to choice 
A in ( |3.4| ) ) . The reduction of the NLL and resummed solutions with respect to the leading 
order solution is substantial in all regions of k. As we see, the resummed solution also has a 
better behavior in the large and small momentum limits as compared to the NLL solution. 
This is due to the absence of the double and triple collinear poles. 

Next, we proceed to the analysis including the saturation boundary. For this, as in 
section pi we use the asymmetric scale choice and the initial condition (|3.1D . The results 



are shown in figure |10[ The resummed solution grows initially rather slow which is to be 
expected as the resummed prescription contains a tower of terms which are subleading for 
high values of rapidity but which are nevertheless important for the phenomenologically 
relevant values of smaller Y. This delay in the growth of the resummed evolution has 
interesting and important consequences for the phenomenology of saturation as it clearly 
implies that the growth of the saturation scale will be delayed in Y. 

That the effects of saturation are thus suppressed in the resummed evolution for lower 
Y is apparent in the right plot of figure ll^ where the nonlinear and linear evolutions 
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Figure 10: Left: The NLL (dotted blue) and resummed (solid red) BFKL evolutions with the 



saturation boundary (2.28). In each case the strong coupling constant is running as as{q) in front 
of the leading order term, and as as(max(fc, k')) in front of the subleading terms. The four sets of 
curves correspond to F = 2, 6, 10, 14. Right: Linear resummed BFKL solution with (solid red) and 



without (dotted blue) the boundary ( 2.28 ), with a running coupling, for the same set of rapidities. 



(i.e. with and without the boundary, respectively) are compared. Clearly, for the values 
of Y shown in the figure, the progress of the front is not much affected by the inclusion 
of saturation. This result can be compared to the earlier results presented for the LL 
and NLL evolutions in figure |7[ It is then clear that saturation has a larger effect on the 
front in the NLL case than in the resummed case. It is interesting that if we look at the 
linear solutions only (the dotted blue lines), then at the highest rapidity in the figure, 
Y = 14, the resummed solution in figure |l^ is actually even slightly larger than the NLL 
solution at the same rapidity in figure ^ Despite this, however, the effect of saturation is 
manifestly smaller in the former. This is due to the fact that the evolution at low rapidities 
is significantly slowed down in the former due to the resummation, and this implies in turn 
that the saturation effects do not become important until later. We also note that the 
resummed evolution is stable and there is therefore no complication or ambiguity in the 
extraction of the saturation scale. From the figures we also observe again the behavior at 
large k which we noticed earlier in the case of the fixed coupling NLL evolution whereby 
the saturated solution actually lies slightly above the linear one. The reason appears to 
be due to the coUinear subtractions which still imply that the solution at very large k can 
turn negative. Unlike the fixed coupling NLL evolution, however, the resummed evolution 
is stable, and its results are unambiguous. Moreover, the resummed prescription is not 
particularly sensitive to the choice of the running coupling, again unlike the NLL evolution 
which displays a high sensitivity. 
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4.2.1 The saturation scale from resummed approach 

It is straightforward to extract the saturation scale from the above results, but let us first 
look at the fixed coupling case where it is easier to make semi-analytic estimates of the 
result. We can then check the consistency of the implementation by comparing these esti- 
mates with the full numerical results. 

Fixed coupling case 

It is possible to calculate analytically the behavior of the saturation scale as a function 
of rapidity in the case of fixed coupling. This was originally done in [33], and later also 
using the traveling wave method in [60]. The resummed evolution was also studied using 
the same methods in [34]. The full expression for fixed coupling is 



Qf(y) = Q^exp -; -'--- y---^-lny . (4.15) 

V 1 - 7s ^(1 - Is) J 

where the "saturation anomalous dimension" 7^ is determined by the condition 

X{ls,Ois) 



1-7. 



X'(7.,a.). (4.16) 



Here the prime on x ™ the right hand side denotes the derivative with respect to 7^. 

Using this formula we can extract the rapidity dependence also for the resummed 
evolution. However, in this case the collinear resummation complicates somewhat the 
calculation because it introduces in the eigenfunction x ^'^ additional dependence on the 
Mellin variable a;, which is itself equal to UsX- That is, we now have a relationship, 

W = Q;sXresum(7,as,w) , (4.17) 

where Xresum was given in ( |4.5| ). This transcendental equation is not analytically solvable 
for uj, but it can be numerically solved. We define a new function XcRil-, c^s) such that for 
any values of 7 and a^, 

UJ = asXea{l,OLs) (4.18) 



gives the value of lo that satisfies equation (4.17). The solution Xeff then fills the role of x 



in equation ( [4.15 ), allowing us to calculate the rapidity dependence of the saturation scale. 
As the absolute normalization of Qs is not under full control it is customary to calculate 
the logarithmic derivative 

It is well known that in the leading order BFKL evolution one has As = 4.88 a^, and that 
js ~ 0.37. 

Constructing the "effective" eigenfunction, Xeflf; numerically, we have calculated A^ for 
different values of Og for the resummed evolution. For rapidities up to 50 units we find 
that Xg is given by 0.308 for ag = 0.1 and 0.528 for as = 0.2. If only the leading term in 



( 4.15 ) is kept then these numbers increase to 0.361 and 0.580 respectively. From the full 



numerical solution we instead find the numbers 0.322 and 0.558 respectively. Notice that 
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Figure 11: The saturation scale obtained using the boundary ( 2.28 ) as a function of rapidity 
for the LL (blue), NLL (green) and resummed (red) evolutions. In each case the strong coupling 
constant is running as as{q) in front of the leading order term, and as as(max(fc, fc')) in front of 
the subleading terms in the resummed and NLL cases. 



these numbers are still rather below the leading order asymptotic result which indicates 
the importance of the non asymptotic corrections even up to relatively large Y. That the 
calculation slightly underestimates the slope obtained from simulation data seems reason- 
able because the highest order corrections that we are ignoring in the calculation are likely 
to be positive. In general, however, the calculation is rather consistent with the simulations. 



The running coupling case and the full numerical solutions 



In figure 11 we present our main result, namely the saturation scale as a function 
of the rapidity for the LL, NLL and resummed prescriptions. We here include only the 
results obtained using the boundary ( |2.28| ) for the resummed case. Let us mention that 
the different boundaries lead to different normalizations of the extracted saturation scales. 
This difference in normalization is due to the evolution only for the lowest values of Y. For 
the leading order solution for example, we find that the solutions obtained from different 
boundary prescriptions have the same Y dependence for Y > 6 (of course the exact numbers 
can depend on the initial condition). For the resummed evolution this number should be 
slightly higher due to the delayed evolution. The overall differences are not that large, 
however. 



We see from figure 11 that the saturation scale obtained from the resummed evolution 
is suppressed at lower rapidities compared to the NLL result, though it should again be 
kept in mind that the NLL solution is generally unstable. In addition, also the absolute 
normalization of Qs cannot be taken too seriously since it depends very much on the exact 
definition of Qs- What is clear, however, is that the resummed result is suppressed for 
the smallest rapidities due to the large-x terms in the evolution. The "plateau" observed 



for the resummed saturation scale at lower rapidities in figure 11 can be attributed to the 
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Figure 12: The logarithmic derivative of the saturation momentum as a function of rapidity for 
the LL (blue), NLL (green), and resummed (red) evolutions obtained using the boundary condition 
( 2.28|) . In each case the strong coupling constant is running as ag (q) in front of the leading order 
term, and as as(max(fc, k')) in front of the subleading terms in the resummed and NLL cases. 



previously mentioned "dip" in the gluon splitting function. While the details of the results 



in figure 11 at the smallest Y inevitably depend on the exact initial condition and also the 
boundary condition (the type of the boundary and the numerical parameters chosen for 
given boundary), we have checked that the plateau of the saturation scale in the resummed 
case also appears in case we change the boundary or when we consider a different initial 
condition (in particular we have considered also an x-dependent initial condition which 
does not appear to wash out the plateau). 

We next show in figure ^ the logarithmic derivative of the saturation momentum, As, 
which was defined in ( [4.19 ). We notice that the RG result is smaller than the LL one in 
the asymptotic region, and that it approaches it slowly from below in agreement with the 
results in [34]. However, for values of Y up to around 15, there is a discrepancy with the 
findings of [34] which appears to be due to the presence of the non-asymptotic contributions 
in our numerical treatment. It is indeed clear from figure ^ that the large-x terms play 
a rather important role for the resummed result, as the red curve makes a sharp turn and 
exceeds the blue curve of the leading order result for smaller rapidities. Even though the 
"speed" of the saturation scale thus defined is relatively large for the resummed evolution 



for Y < 10, notice that the saturation scale itself is small as is evident from figure 11. 

We thus see that the preasymptotic behavior of the resummed solution has important 
consequences for the study of saturation. How far exactly the plateau extends in Y is 
hard to answer in our present application of the saturation boundary. To answer this 
question exactly, we would need to perform a serious phenomenological study where all 
the unknown parameters would be fitted to data. Regardless of the exact values of the 
parameters, however, it seems to be clear that resummation is needed for sensible results^. 



^It is, strictly speaking, not a problem if the gluon Green's function goes negative at some k since 
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and that it implies that the growth of the saturation scale is delayed. 

5. Summary 

In this work we have studied small-x evolution in the presence of saturation beyond the 
LL order. We have performed our analysis by self-consistently solving the linear BFKL 
evolution with proper boundary conditions; this is a justified procedure so long as we are 
interested in determining the energy dependence of the saturation momentum Qg, and 
any observable under consideration for momenta above (but also close to) Qs- We have 
found that the NLL corrections are large and that it is therefore not enough to keep 
only those next-to- leading terms related to the running of the coupling. Moreover, and as 
expected from a direct inspection of the NLL BK equation or from the well-known behavior 
of the NLL BFKL equation, the corrections can become negative in certain regimes of 
the transverse (momentum or coordinate) space and the NLL analysis is unstable. Thus 
we have used the renormalization group improved evolution kernel which resums all the 
dominant N^LL contributions arising in the collinear (and anticollinear) limit, for all n > 2, 
and has a stable behavior. In accordance with previous expectations we do find that at high 
rapidities the running coupling results for the RG improved evolution converge to those 
of the LL evolution for the logarithmic derivative of the saturation momentum. However, 
for smaller values of the rapidity, that is for values up to y = 10 or 15, we find significant 
deviations from the asymptotic and the pure running coupling results. At least for the 
initial conditions studied, we find that the saturation momentum does not grow in the first 
few units of the rapidity interval, and this is presumably the manifestation in the context 
of saturation of the well-known "dip" of the gluon-gluon splitting function. From that 
point on, Qg increases rather fast for the next few units until it becomes reliable to use 
the asymptotic expressions. At this stage, we do not know if our result is robust for all 
types of hadrons, e.g. for a large nucleus, and if it is likely to change when one properly 
includes the relevant impact factors. However, it is an interesting phenomenon which, to 
our knowledge, has not been studied so far and clearly deserves further investigation. 
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